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ABSTRACT 

A mathematical model has been developed for spinning mode acoustic 
radiation from a thick wall duct without flow- This model is based on a 
series of experiments (with and without flow) conducted by Richard Silcox 
[Silcox] of the Noise Control Branch at Langley Research Center. In these 
experiments a nearly pure azimuthal spinning mode was isolated and then 
reflection coefficients and far field pressure (amplitude and phase) was 
measured. In our model the governing boundary value problem for the Helmholtz 
equation is first converted into an integral equation for the unknown acoustic 
pressure over a disk, SI, near the mouth of the duct and over the exterior 
surface, S2, of the duct. Assuming a pure azimuthal mode excitation, the 
azimuthal dependence is integrated out which yields an integral equation over 
the generator Cl of SI and the generator C2 of S2 (see Figure 2) . We 
approximate the sound pressure on Cl by a truncated modal expansion of the 
interior acoustic pressure. We use piecewise linear spline approximation on 
C2. We collocate at the knots of the spline and at zeros of the first term 
excluded in the truncated modal expansion. Finally, we compare numerical and 
experimental results . 


Research was supported by the National Aeronautics and Space Administration 
under NASA Contracts No. NAS 1-1 7130 and NAS 1-1 6394 while the author was in 
residence at the Institute for Computer Applications in Science and 
Engineering, NASA Langley Research Center, Hampton, VA 23665. 
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INTRODUCTION 


The reduction of jet noise radiated from the inlets of aircraft engines 
has become an increasing important problem. Methods to suppress this noise 
have included the development of acoustic liners, high Mach number inlets, and 
the use of inlet geometry to redirect the sound. Experiments with and without 
flow have been conducted at Langley Research Center [Ville] , [Silcox] to study 
these methods. In particular, in July 1982 an experimental study was 
conducted by Richard Silcox [Silcox] in a continuing effort to examine the 
effect of inlet geometry on the reflected and radiated acoustic fields. 

This report describes a mathematical model for the no-flow experiments. 
This model has been constructed to cover a range of reduced wave numbers and 
azimuthal and radial modes. The inlet contour can be modified so that the 
result is a model for both inlet and engine, each with its own excitation. An 
axisymmetric impedance boundary condition can also be easily incorporated into 
a portion of the interior duct wall. 

The experiments of Richard Silcox were designed around the spinning mode 
synthesizer (SMS) in the Langley Research Center flow duct facility. The SMS 
can excite a nearly pure (20-30 dB isolation) azimuthal mode inside the 
duct. This modal expansion for the interior acoustic pressure is built into 
the model. At an artificial circular disk interface, SI, near the mouth of 
the duct, the pressure and its normal derivative are required to match 
continuously (see Figure 2). On the exterior duct surface, S2, a hard wall 
boundary condition is imposed. In the exterior region, the pressure is 
required to satisfy the Helmholtz equation and the radiation condition at 
infinity. This boundary value problem is converted into an integral equation 
over SI + S2 using Helmholtz' formula. The unknowns in this equation are 
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the complex pressure on S2 and the reflection coefficients in the interior 
modal expansion. 

By assuming a single azimuthal mode excitation, it is possible to 
integrate out the azimuthal dependence which yields an integral equation over 
the generator Cl + C2 of SI + S2. We arrive at a one-dimensional integral 
equation at the expense of a somewhat more complicated kernel. 

The numerical method used is collocation; i.e., the integral equation is 
required to be satisfied at certain points on Cl 4- C2. The unknown pressure 
on Cl (and its normal derivative) is approximated by a finite Bessel series 
which is a truncation of the interior azimuthal mode expansion. Piecewise 
linear spline approximation is used on C2. The absolute error in the 
solution is estimated at the knots of the spline and this information is used 
to recommend the number of knots required for a given error tolerance on 
C2. This information can also be used to distribute the recommended number of 
knots to achieve an equal distribution of the absolute error. This is useful 
if it should be necessary to run the code a second time in order to achieve a 
better error performance. The code also provides (optionally) for one step of 
Neumann iteration. This yields a natural interpolation formula for the 
pressure, and gives an approximation with the same smoothness as the exact 
solution. Finally, Helmholtz' formula is used to compute the pressure on a 
semicircle in front of the duct for comparison with experimental results. 

In Section 2 we give for completeness a brief description of the SMS and 
the Langley Research Center flow duct facility. In Section 3 we describe the 
mathematical details of the model, while in Section 4 we outline the numerical 
method. Our codes and their implementation are discussed in Section 5. 
Numerical results are presented in Section 6 and Section 7 contains concluding 


remarks . 
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2. THE SPINNING MODE SYNTHESIZER (SMS) 

For completeness, we give a brief description of the SMS. The reader is 
referred to Figure 1 and references [Pal] and [Ville] . The following is taken 
from [Ville] . 

The SMS incorporated into the flow duct facility is a research apparatus 
designed to overcome the problems involved in static testing of real turbofan 
engines or research fans. The SMS generates arbitrary combinations of 
acoustic patterns at a specified frequency in a 0.3 meter duct (with or 
without airflow). Specified duct modes are generated by controlling the 
amplitude and phase of 24 acoustic drivers equispaced around the duct wall in 
a plane perpendicular to the duct center line* By properly adjusting the 
input to the drivers, individual spinning modes, a combination of modes, or 
circumferential standing waves can be generated in the duct. The acoustic 
field produced by the array of drivers is monitored by an array of 48 wall 

mounted microphones located 0.2 meters upstream of the drivers. At these 

microphone locations the desired acoustic wall pressures (amplitude and phase) 
are approximated to some specified degree of accuracy. In order to attain 
this accuracy, the pressure field sensed by the 48 microphones feeds back 

through a control computer optimization algorithm to generate correction 

signals to the drivers. By an iterative process, the pressure field at the 
microphone array converges to some specified target pressure at the monitoring 
array. Once the target pressure is attained, the mode setup may be stored by 
recording the driver settings for future recall, or the steady state acoustic 
field may be left intact for experimental studies. 

The inlet duct test apparatus is mounted in the flow duct facility of the 
Langley Aircraft Noise Reduction Laboratory and extends into the anechoic 
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chamber shown by the view in Figure 1. This chamber measures 9.15 meters 
wide, by 6.1 meters deep by 7.16 meters high. 

The flight inlet is attached to the inlet section of the duct. The 
geometry actually used in our model is indicated in Figure 2. We have 
extended the flight inlet contour BCD with a straight (r = constant) portion 
DE, and quarter arc of a circle EF. We note that our numerical results 
indicate that this termination of the inlet has little effect at the 
frequencies considered . 


3. THE MODEL 

This model was suggested by the approach of Kagawa, et al. [Kag] to 
loudspeaker design. 

Figure 2 gives a contour drawing of the idealized inlet. This is the 
flight inlet contour (an arc BC of an ellipse plus the line segment CD) 
with a straight extension DE and a quarter circle termination EF. 

Let a denote the interior duct radius (a « 0.15 meters) and introduce a 
cylindrical coordinate system (z,r,0) with origin 0 and positive z axis 
pointing out of the duct. Let m denote angular frequency, c the speed of 
sound, and k = w/c the reduced wave number. We will use the dimensionless 
coordinates z = z/a and r = r/a and the parameter k = ka. Let ft 
denote the exterior of SI + S2 and let D denote the interior of SI + S2. 

We denote complex pressure by $(z,r,0)e where 

($, -d < z < 0, r < 1 

*- 

U 


in ft 
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For a fixed integer m, it is assumed that 

— im6 , u 
expansion <p = <pe with 


the SMS generates the modal 


(3.1) <f> = l (A(n)e lL ^ Z + R(n)e iL ^ ) J m (X(n)r). 

n=0 

This expansion comes from separation of variables in the reduced wave equation 

plus application of the hard wall boundary condition. In practice other m 

or azimuthal modes are present inside the duct but are at least 20 dB below 

the desired mode. In (3.1) J ffi denotes the ordinary Bessel function of the 

first kind and order m. The increasing sequence X(n) is defined by 

J'fX(n)l = 0, n > 0. We assume that k * X(n) for all n and define NCT 
nr J 

to be the integer satisfying X(NCT-l) < k < X(NCT). Then 


L(n) 


= / k 2 - (X(n)) 2 , 


0 < n < NCT-1 


(3.2) 


-iL(n) 



n > NCT. 


The radial modes corresponding to n=0, • • • , NCT-1 are called cuton, the other 
cutoff, and NCT is the number of cuton radial modes. Complex R(n) is 
called the reflection coefficient of the n th radial mode, while complex 
A(n) is the strength of the forward propagating modes, since we assume that 
A(n) = 0 for n > NCT. This is a reasonable assumption since the plane of 
the 24 acoustic drivers of the SMS is about 13 duct interior diameters from 

SI. 

Note that the dependence on m in the above notation has been suppressed. 
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The appropriate boundary value problem for the Helmholtz equation (the 
wave equation with the harmonic time dependence separated out) can be stated 
as follows: 

Find ^ in class C^(fi) c\fi) such that 
Aij; + k 2 t = 0 in ft, 

3i|i 

(3.3) -gj = ° on S2, 

t satisfies the radiation condition at infinity, 

* = ♦ and = on S1 * 

We use Helmholtz' formula (an application of Green's third identity) to 
convert (3.3) into an integral equation. Let £ denote the observation point 
and the integration point. Then 



where R = lx - x' I • Here tT denotes the normal to SI + S2 pointing 
into D. Now uniqueness for (3.3) implies that = ^(z,r)e^ in ^. If we 
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change to cylindrical coordinates in (3.4), multiply by e , and integrate 
from 0 = — 7T to 9 = tt, we obtain the one-dimensional integral equation 


(3.5) 


^(z,r), (z,r) € u - (Cl n C2) 

?(z,r)/2, (z,r) = Cl (1 C2 

<t>(z,r) , (z,r) € Cl 


= - / r' dr' f-r <K0,r') K(z,r;0,r') 

o 8z 


+ / r' dr' <K0,r') -gpr (z,r;0,r') 


- I 

C2 


3 V 

r' ds' t(z',r') -g^r (z,r;z',r')» 


where 


(3.6) 


it ikR 

K(z,r;z',r') =ir / cos(mY) ~ dY, 


0 


R 


(3.7) 


R 


= ((z-z') 2 + (r-r') 2 + 4rr' sin 2 (Y/2) ) ^ 2 , 


s' denotes arc length on C2, and n' denotes the normal to C2 pointing 
toward the duct interior. We note that (3*5) appears to be homogeneous; 
however, <{> in (3.1) is the sum of two terms, one of which is assumed known. 
This is the excitation term 


| A(„). iI - <n)2 J (i(n)r) 

n=0 


controlled by the SMS. 
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4. THE NUMERICAL METHOD 

Our numerical method is collocation. We require (3.5) to hold at 
specified (collocation) points on Cl + C2. We obtain a full, square complex 
linear system whose solution provides an approximate solution to (3.5). 

For the approximation of <|> we truncate (3.1) at n = N1 > NCT-1, 

N1 

(4.1) 4> - t (A(n)e iL ^ n ^ Z + R(n)e -iL ^ n ^ Z ) J (X(n)r). 

n=0 m 

To approximate t on C2, we first parameterize this contour. We will refer 
to certain values of this parameter t as knots. We require points B, C, D, 
E, and F to be knots. Additionally, knots are added so that BC is divided 
into K2 sub intervals of equal length with respect to this parameter, CD 

into K3 subintervals, DE into K4 subintervals, and EF into K5 

subintervals. This yields N2 + 1 = K2 + K3 + K4 + K5 + 1 knots. 

An additional line segment AB can be added to the duct mouth and divided 
into K1 sub intervals. This is useful if it is desired to alter the hardwall 
condition on AB. In the following A = B and K1 = 4*, but activation of 

these parameters is provided for in our codes. 

We use N2 Chapeau functions, t n (s) depending on arc length s to 

approximate These basis functions are centered at the values of s 

corresponding to the knots, except for the knot corresponding to the endpoint 
F of C2. This is because the solution at F must be zero for azimuthal 
mode index m > 1. Hence, we write 

_ N1+N2 _ 

(4.2) ip = l C(n) i|> (s) on C2. 

n=Nl+l n 
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We have experimented with two numerical procedures. 

I. (i) Collocate at (0,r Q ) , • •* , (0,r N1 ) where 0 < r Q < ••• < r N1 < 1 

satisfy J m ( A(Nl+l)r^) = 0; i.e., choose the r coordinates to be 
positive zeros of the first term left out in the truncation (4.1) 
with z = 0. 

(ii) Collocate at the points ^ Z N1+1 * r Nl+l ^ * * * * * ^ z Nl+N2 ,r Nl+N2^ 

corresponding to the centers of the Chapeau functions. 

II. (i) as in I 

(ii) as in I but replace the equation generated by the collocation point 
at B by the continuity equation 

N1 N1 

(4.3) - l R(n) J (X(n)) + C(N1+1) = l A(n) J U(n)). 

n=0 m n=0 

Now let R = [R(0), •••,R(N1)] T , C = [ (Nl+1) , • • • ,C(N1+N2) ] T , and 
A = [A(0) , • • • ,A(N1) ]^. We write the linear system resulting from I or II by 

(4.4) CKT0T*[R,C] T = CRITE*A , 

where CKTOT is a full, complex square matrix of order N1 + N2 + 1 and 

CRITE is an (N1 + 1) * (N1 + N2 + 1) complex matrix. Writing (4.4) in 

block form for procedure I we have 

(4.5) 


] ["&] p-j*n 
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where 


(4.6) = J m (X(n)r^), 0 < £,n < Nl, 


(4.7 


, 0 C £, n < Nl 

S£(l,n) y Nl+1 < £ < Nl + N2, 0 < n < Nl 

1 


-iL(n) ^ r' dr' J m (X( n )r') k(z £ . r^ ;0,r') 


(4.8) i*U,n) = -J r' dr' J m (X(n)r') ffr (z r r A ;0,r'), 


Nl+1 < £ < N1+N2, 0 < n < Nl, 


(4.9) 


#(£,n), 0 < £ < Nl, Nl+1 < n < NI+N2) 

4f(£,n), Nl+1 < £,n < N1+N2 ) 


- / r' ds' t (s') 

C2 n 3n 


- (z. ,r. ;z' ,r') , 


^(Nl+1, Nl+1) = I/2 

(4.10) <^(£,£) 1, Nl+2 < £ < N1+N2 

^(i,j) = 0, Nl+1 < i*j < N1+N2. 

We note that F(£,n) = 0 for 0 < £ < Nl+1 because — 7 (0,r:0,r') = 0 for 

o Z 

0 < r * r' < 1. 
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If procedure II is used, then row N1 + 1 in CKTOT and CRITE must be 
changed according to the continuity equation (4.3). 

The system (4.4) is solved NCT times with "basis" excitations 
Aq, • • • >A N ct_i which satisfy A^(i) = 6^. We cal] the resulting approximate 
solutions to (3.5), basis approximate solutions and denote them by end 

i 1 * 

Calculation of the basis far field (or basic first Neumann iterate) is 
done using the Helmholtz' formula: 


(4.11) 


2, (z,r,0) in 

F"(z,r) •{ 1, (z,r) € Cl + C2 - (Cl fl C2) 

V 2 » (z>r) = Cl 0 C2 


1 3j> 

- / r' dr' v^(0,r') K(z,r;0,r') 
0 9z 


+ / r' dr' /(0,r') ^ (z,r;0,r') 
0 


9K 

3z" 


-Jr' ds' ^(z'.r') (z,r;z',r'). 

C2 1 


Given the coefficients A(0) , • • • ,A(NCT -1) of the desired excitation, we find 
the corresponding approximate solution to (3.5) and the far field (or first 
Neumann iterate) from 


(4.12) 


NCT— 1 NCT-1 

i- I A(iH\ F = l A(i)F\ 
i=0 i=0 


The initial choice of N1 and N2 in (4.1) and (4.2) may not always be 
consistent with the desired accuracy. Also, the equispacing of the knots with 
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respect to the parameter t is usually suboptimal. Section 6 discusses rules 
for choosing N1 and N2 initially. We now present a procedure (see 
[deBoor2]) which uses the first approximate solution obtained to estimate the 
absolute error at the knots, and then recommends new values for K2, K3, K4 
and K5 and a new distribution of the knots. The goal being to equally 

distribute the absolute error among the knots and to achieve a desired 
accuracy. 

Let us explain this procedure for BC which is initially divided into 

K2 sub intervals . Let t^ denote an interior knot and let h^ and h i+1 

denote the mesh spacing immediately to the left and right of t^. Set 

h = max{h^,h^ + ^}- At t^ we estimate the second derivative with respect to 
th 

t of the j basis approximate solution to (3.5) by interpolation with a 
parabola at t^_p t^ and t i+l* Denote this estimate by d^ and set 
d^ = max(|d^|: 0 < j < NCT-l}. We estimate the absolute error at t^ by 
d^ h^/2. Let AE denote the sum of the errors at the interior knots divided 

by the number of interior knots, K2— 1* Let AER denote the desired absolute 
error, and set 

L2 = [K2//AER/AE] + 1, 

unless the result is 1, in which case set L2 = 2. 

Next, we partition BC into L2 sub intervals so that the error is 
approximately the same at all interior knots. Set p = K2 and l = L2 and 
let denote the initial knots on BC. Set u i “ u ± ~ 

(t^ + i = 2,* ## ,p-l, ^ = tp + ^. Let G denote the piecewise 

constant function defined by G(t) = /d^ + ^ for u^ < t < u i+ j. For 
u^ < u < Up, define 
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u 

u p 

y = / G(t)dt and y = / G(t)dt. 


u n 


u 


1 


Now y is increasing in u so let u = H(y) denote the inverse function. 
Construct £ + 1 new knots according to = tj, = h( (i-1 , 

i = 2, •••,£, = t p+i* This procedure is based on our observation that in 

this problem the redistribution of a fixed number of knots has little effect 
on the average error AE. 


5. CODE DESCRIPTION 

We will discuss the following codes: CLLEQS (CLLIQS), CLLRDS (CLLIRS), 

NELMAN, FALFLD, and SPDATA. CLLEQS use the continuity condition (4.3) at 
Cl C2, while CLLIQS and CLLIRS collocate at this point. Our codes are 
programmed in CDC Fortran 5 (Fortran 77). 

We begin by discussing CLLIQS. First, values for k = ka (RKA), m, Nl, 

N2, NCT, and AER, the error tolerance for C2, must be input. CLLIQS contains 

the data X(0) , • • • , X(NIMAX) , X(0) , • • • , X (N1MAX) in arrays ZERO and RLAMDA for 
m = 1,***,4 (currently), where J m (X(i)) = 0 and J^(X(i)) = 0. Array RL 
contains the complex parameters L(0) , • • • ,L (N1MAX) . Currently, N1MAX = 13. 
The value of NCT, the number of cuton modes is computed from RLAMDA and k, 
and then checked against the user supplied value. 

Subroutine CURV contains the parameterization of C2 and provides, for a 

given parameter value t, the corresponding value of z and r, the z and 

r components of the normal, and the value of arc length measured from B. 

The values of the constants determining the location of B, C, D, E and F are 

stored in subroutine CONSTS. 
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BC is an arc of the ellipse in the z— r plane given by 

f z ->2 rr_^_J^ 1 ^807JL9'i 2 _ . 

.4161438' *• .2080719 J ~ i ’ 

which joins B => (0,1) to C = (.147804, 1.4025775). CD is the line segment 
joining C to D = (-1.6187892, 1.7381843). DE is the line segment 
® bo E = (—2.6, 1.7381843), and EF is an arc of a circle of 

radius 1.7381843 with center (-2.6, 0) joining E and F = (-4.3381843, 0). 
Arc length on BC is computed by approximating 

t 

S(t) = .4161438 / /(5 + 3 cos(2t))/ 8 dx 

0 

using 10 point, piecewise cubic Hermite interpolation for 0 < t < 2.7784911. 
This is accomplished in subroutine HERM and BARCL. 

Subroutine KNOTGN generates the knots (array RKNOT) on C2 in terms of 
the parameter t as described in Section 4. Subroutine COLLGN uses CURV and 
KNOTGN to produce the N1+N2+1 collocation points (arrays COLLZ, COLLR ) and 
the values of arc length corresponding to the knots (array ARCL) . 

Subroutine COEFF assembles the complex matrices CKTOT and CRITE in (4.4). 
Almost all the execution time of CLLIQS is devoted to this one subroutine. 
The difficulty is that the oscillatory and singular double integrals in (4.7) 
- (4.9) must be computed. We use adaptive integration which is ideal for this 
type of behaviour, but expensive. We were constrained by the need to test our 
formulation and produce a working code during 14 weeks at ICASE in the summer 
of 1982. Alternatively, we could have developed a suitable product 
integration formula and then used the Nystrom method instead of collocation, 
but this would have required more time. 
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We have separated the integrands in (4.7) - (4.9) into a bounded part and 

a singular part corresponding to the axisymmetric potential equation. 

Let p = ((r'-r ^) 2 + (z'-z ^) 2 ) ly ^ 2 and P = ((r'+r^) 2 ^- (z'-z , let 

C = kR/2, n' = (N z ^, N^,), and let % and j# 7 denote the complete elliptic 

integrals of the first and second kinds as functions of the complementary 

2,-2 

parameter m^ = p /p . 

Then in (4.7) we have 


^(A,n) = - 


iL (n) 


1 


^ / r' dr' J (X(n)r') / dY(cos(mY) (e ikR -l) + cos (mY)-l)/R 
0 m 0 


_ iL(nI j r , dr , j (X(n)r') / 
17 n m n 


dY 


I TT 

= - — / r' dr' J m (X(n)r') / dY(i cos(mY) sin(?)e i ^ - sin^(mY/2)) /£ 


2iL(n) 


- — / r' dr' J m (X( n )r') J^On^/p , 


while in (4.9) we have 


n) = 


*jy / ikR 

tt - 1 / r' ds' t (s') / dY [cos (mY) 7 (- — + (cos(mY)-l) 7 £ 
C2 n 0 \ an R an R 


*) 


+ it " 1 / r' ds' * (s') J dY f -7 £ 
C2 n 0 an R 


a 1 


= ^ J r' ds' + n (s') J dY^N z ,(z'-z £ )+N r .(r'- rjl + 2r & sin 2 (Y/2))j 

i cos(mX)e U k C ° S ; 2 ~ sln C + i + | 3l °<f /2) j j 
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+ IT -1 / ds' ~ 
C2 p 


^( 8 ') 


2r'jf(m 1 )(N r ,(r i -r') + (z £ -z' ) ) /p 2 + J^))! 


For the evaluation of Jj, J m , m > 2, and we use the LRC 

FTN5MLB library routines BJ1R, BKIR, ELIPKC and ELIPEC. For the evaluation of 
the one-dimensional integrals we use the FTN5MLB routine CADRE which is a 
modification of an algorithm due to deBoor [deBoorl] . For the double 

integrals we use the FTN5MLB routine CAREDB which computes the integral as 
iterated single integrals with the single integrals computed as in CADRE. 

CADRE is an adaptive cautious Romberg extrapolation routine which is 

designed to identity certain types of integrand behaviour by examining a ratio 

based on the previous three trapezoidal sums . We split the integrals so that 

the singularities are endpoint singularities. If such a singularity is 

2 

detected, CADRE switches to a process similar to Aitken's 6 process to 

estimate the integral and evaluate the error. As a result of this switching, 
we have observed that execution times decrease with the error tolerances. We 
use EPS(l) = EPS (2) = E-5, where EPS(l) is the maximum allowable relative 
error and EPS (2) the maximum allowable absolute error. EPS (3) is the sum of 

the computed error estimates over each subinterval in which the convergence 

criteria is satisfied. These remarks also hold for CAREDB. 

Adaptive integration is particularly suited to oscillatory and (or) 
singular integrals. It offers the possibility of constructing a test code for 
an integral equation method in a relatively short time. Both CADRE and CAREDB 
require real integrands so real and imaginary parts were integrated 
separately. We did not take the time to convert these routines into 
integrators for complex valued functions. 
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Once the linear system (4*4) is assembled, we solve it using FTN5MLB 
routine CXGCOIT. This subroutine performs an LU decomposition, solves by 
forward and backward substitution, and estimates the condition number 
(CONDNUM) of CKTOT in the 1-norm. Optionally, iterative refinement can be 

performed (see [Don] ) . 

Finally, recommended values for K2, K3, K4 and K5, with respect to an 
absolute error tolerance (AER) on C2 are computed, and estimates made of the 
second derivative of the approximate solution (with respect to the parameter 
t) at the interior knots of BC, CD, DE and EF. 

On execution a listing is provided of all input parameters, NCT, RLAMDA, 
ZERO, COLLZ, COLLR, RKNOT, ARCL, CONDNUM, the real and imaginary parts of the 
complex approximate solution (as well as phase and amplitude) , a continuity 
check at Cl fl C2, recommended values for K2, K3, K4, K5, and estimates of 
second derivatives and error on C2. In addition, the maximum integration 
error estimate and the corresponding indices is written for the matrices 

and for J^, and tor and 

Program CLLRDS is essentially the same as CLLEQS except that provision is 
made in subroutine KNOTGN for using the recommended values of K2, K3, K4, K5 
and the second derivative estimates to distribute this new number of knots. 

Program NELMAN takes the output of CLLEQS, CLLIQS, CLLRDS or CLLIRS and 
computes one basis Neumann iterate. This iterate and corresponding 
integration error estimates are output. 

Program FALFLD takes the output of CLLEQS, CLLIQS, CLLRDS, CLLIRS or 
NELMAN and computes an approximate basis far field. The basis far field and 
corresponding integration error estimates are output. 
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Program SPDATA takes the output of FALFLD and user provided excitation 
data A(0) , • • • ,A(NCT-1 ) (amplitude in dB, phase in degrees) and computes and 

outputs the approximate solution and far field. The far field amplitude is 
referenced to its maximum value. Finally, plotting files are created. 


6. NUMERICAL RESULTS 

All computations were performed on a CDC Cyber 175 at NASA Langley 
Research Center. We give results for the programs which collocate at Cl fl C2 
instead of demanding the continuity equation (4.3) hold. We have found little 
difference in the outputs of these two sets of programs. 

For a given azimuthal mode index m, it is best to start with a value of 
k = ka so that one mode is cuton (NCT = 1) and then increase k to the 

desired level. Good starting values then are N1 = NCT + 5 and K2 = 20, 

K3 = 10, K4 = 10, K5 = 7. Program CLLIQS recommends new values for K2,***,K5 

with respect to a user supplied tolerance. The user has the option of running 

CLLIRS with these new values and derivative estimates from CLLIQS as input. 
We have provided a continuity check which outputs the modulus of the 
difference in the two sides of (4.3) for an approximate solution. If this 
continuity check is greater than the error tolerance, then N1 can be 

increased in CLLIRS. 

The error estimates on C2 and the continuity check are a reasonable 

indication of the accuracy attained. Of course, the condition number and the 
accuracy with which CKTOT and CRITE are computed must also be considered. We 
have run CLLIRS for k = ka = 2.66 and an increasing sequence of values for 
N1 + N2. Table I gives some of these results. CONT is the value of the 


continuity check. The error in |R (0) | and the error at Cl fl C2 are 
estimated by comparison with the case N1 + N2 = 78. These data support the 
hypothesis that the condition number in the 1-norm satisfies 
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CONDNUM < 1.8 N1 + 4.5 / N2 , 


and that the error in R(0) and the maximum error on C2 are less than 


max 




for some constants and I^. 

Table II gives for m = 1 a comparison of experimental and numerical 
results. Some of these results are in close agreement; for example, for 
k = 2.66, we have |R(0) |/|A(0) | = .348 (numerical) and .35 (experimental). 
However, for k = 3.20, we have |R(0) | / 1 A(0) | = .164, while the experimental 
value is .196. Our error tolerance here is .002. The continuity check and 
error estimates on C2 are consistent with this tolerance. Generally, when 
R(i) is small compared to max{ | A(0) | , • • • , | A(NCT-l) | } , we have the greatest 
relative error. We have indicated the error bracket on some of these entries 
in Table II. But these brackets do not account for some of the discrepancies. 
We note that the dB levels in Table II are referenced so that the peak sound 
pressure level in the far field is 0 dB (for both experimental and numerical 
results) . 

Figures 3, 4, 5, 6, and 7 give far field patterns computed on a semicircle 
of radius 20 in front of the duct (see Figure 1). We note the good agreement 
between the numerical and experimental curves. The small oscillations in the 
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experimental curve may be due to reflections, while the lack of symmetry 
indicates less than complete isolation of the desired azimuthal mode. We 
compute and store a basis approximate solution and far field. Then we can 
interactively produce approximate solutions and far field patterns for any 
given excitation strengths A(0) , • • • , A(NCT-l) • Thus we have an NCT parameter 
family of possible far fields that can be generated quickly. It is possible, 
for example to make the -30 dB downward spikes in the k = 6.50 case almost 
completely disappear by choice of excitation strengths. But roughly speaking, 
these different patterns have the same "envelope." 


7. REMARKS 

We will discuss existence, uniqueness and regularization in detail in a 

later report. For now we note that the standard proof using a theorem of 

Rellich (see, i.e., [Hellwig] ) shows that uniqueness holds for the boundary 

value problem (3.3). On the other hand uniqueness may not hold for the 

integral equation (3.5). If uniqueness does not hold for (3.5), it follows 

that k is an eigenvalue for Helmholtz' equation in the interior of SI + S2 

with Dirichlet boundary conditions, and that the corresponding eigenfunction 

im 0 

has the form u(z,r)e . If k is such an eigenvalue, we have not been 
able to adapt the usual method to produce a homogeneous solution to (3.5). 
However, the numerical evidence, namely the variation of conditioning with 
changes in k suggests that nonuniqueness does hold at such interior 
eigenvalues. 

We have considered the problem of regularizing (3.5). A good discussion 
of this problem can be found in [Burton] and the references there. We have 
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developed a regularization procedure based on a vector relationship due to 
Maue (see [Burton]). This procedure uses the Galerkin method, Maue's formula 
and Stokes' theorem to move a derivative from the kernel to the unknown 

function. This reduces the strength of the singularity which must be 

integrated. This procedure has not yet been fully implemented and tested. 
Our current code provides condition number estimates, a continuity check (when 
collocation is used at Cl fl C2), and error estimates on C2. These should 
indicate when the results are unreliable. 

We have experimented with the best location for the interface surface 

SI. This is a trade-off between the two types of approximation - Bessel 
series and piecewise linear splines. The best efficiency is obtained by using 
the modal expansion over as large a region as possible. 

We have already mentioned the idea of a full engine model* This more 

complex geometry would be stored in subroutine CURV. Two interface surfaces 
and two modal expansions would be required. This is a natural extension of 
the current code. It is also possible to install a variable geometry. The 
input would be points on the outer duct contour and the code would 
automatically connect these points with a smooth curve having some specified 
property. 

In conclusion we point out that duct acoustics has been studied 
extensively in recent years. The reader is referred to [Baum] for a 
bibliography and a discussion of the various numerical methods which have been 


used. 
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Table I 


N1 

N2 

CONDUM 

CONT 

|R (0) | 

ANGLE 

R(0) 

ERROR 

|R(0)| 

ERROR 
Cl 0 02 

MAX 
ERROR 
EST. C2 

6 

18 

29 

.0016 

.3555 

-20.76 

.011 

.0077 

.0080 

2 

34 

26.8 

.006 

.3479 

-20. 13 

.0029 

.0012 

.0021 

6 

34 

34.1 

.00076 

.3476 

-20.17 

.0027 

.0015 

.0019 

10 

34 

41.2 

.00044 

.3475 

-20.17 

.0026 

.0019 

.0019 

6 

68 

47.9 

.0011 

.3451 

-20.00 

.000043 

.00037 

.00051 

10 

68 

54.2 

.00038 

.3451 

-20.00 





Table II- Max Far Field SPL = 0 dB 


NUMERICAL 

EXPERIMENTAL 

k = ka 

Coefficient 

dB 

Angle 

dB 


A(0) 

36.0 

0° 

36.1 

2.66 

R(0) 

26.8 

O 

CM 

• 

o 

CM 

1 

27.0 


A(0) 

33.8 

0° 

32.6 

3.20 

R(0) 

18.1 ± .2 

25.2° 

18.5 


A(0) 

24.4 

-8.4° 

25.8 


R(0) 

2.2 + 1.64 

29.8° 

18.3 

5.54 

A(l) 

37.8 

-139.1° 

39.2 


R(l) 

28.5 

141.3° 

29.6 


A(0) 

23.9 

70.4° 

22.9 


R(0) 

-9.2 ± 4.3 

19.90° 

9.9 

6.50 

A(1 ) 

34.5 

-64° 

33.5 


R(l) 

-1.3 ± 1.7 

-54° 

9.9 


A(0) 

23.4 

-159.6° 

23.5 


R(0) 

-5.4 ± 2.2 

-168.2° 

3.4 

7.68 

A(l) 

33.0 

56.6° 

33.1 


R ( 1 ) 

-4 ± 1.2 

-5.7° 

9.5 
















































1. Report No. 

NASA CR-172273 

4. Title and Subtitle 


2. Government Accession No. 


3. Recipient's Catalog No. 


5. Report Date 
November 1983 

Spinning Mode Acoustic Radiation From The Flight Inlet I — g n , : — : — : — - — 

° 6. Performing Organization Code 


15. Supplementary Notes 

Langley Technical Monitor: Robert H. Tolson 

Final Report 

16. Abstract 

A mathematical model has been developed for spinning mode acoustic radiation from 
a thick wall duct without flow. This model is based on a series of experiments 
(with and without flow) conducted by Richard Silcox [Silcox] of the Noise Control 
Branch at Langley Research Center. In these experiments a nearly pure azimuthal 
spinning mode was isolated and then reflection coefficients and far field pressure 
(amplitude and phase) was measured. In our model the governing boundary value 
problem for the Helmholtz equation is first converted into an integral equation for th< 
unknown acoustic pressure over a disk, SI, near the mouth of the duct and over the 
exterior surface, S2, of the duct. Assuming a pure azimuthal mode excitation, the 
azimuthal dependence is integrated out which yields an integral equation over the 
generator Cl of SI and the generator C2 of S2 (see Figure 2) . We approximate the 
sound pressure on Cl by a truncated modal expansion of the interior acoustic pressure. 
We use piecewise linear spline approximation on C2. We collocate at the knots of the 
spline and at zeros of the first term excluded in the truncated modal expansion. 
Finally, we compare numerical and experimental results. 




For sale by the National Technical Information Service, Springfield, Virginia 22161 

NASA-Langley , 1984 


8. Performing Organization Report No. 

83-63 

10. Work Unit No. 

11. Contract or Grant No. 

NAS 1-16394 , NAS1-17130 

13. Type of Report and Period Covered 

contractor report 

14. Sponsoring Agency Code 


7. Author(s) 

William F. Moss 

9. Performing Organization Name and Address 

Institute for Computer Applications in Science 
and Engineering 

Mail Stop 132C, NASA Langley Research Center 
Hampton, VA 23665 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546 

















